02_BayesianLinearRegression (Score: 5.25 / 5.0)

  1. Coding free-response (Score: 2.0 / 2.0)
  2. Coding free-response (Score: 1.0 / 2.0)
  3. Coding free-response (Score: 1.0 / 0.0)
  4. Coding free-response (Score: 1.25 / 1.0)
  5. Comment

Przed oddaniem zadania upewnij się, że wszystko działa poprawnie. Uruchom ponownie kernel (z paska menu: Kernel$\rightarrow$Restart) a następnie wykonaj wszystkie komórki (z paska menu: Cell$\rightarrow$Run All).

Upewnij się, że wypełniłeś wszystkie pola TU WPISZ KOD lub TU WPISZ ODPOWIEDŹ, oraz że podałeś swoje imię i nazwisko poniżej:

In [1]:
NAME = "MICHAL MARCINIAK"

In [2]:
import pandas as pd

from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn import linear_model

import pickle
import numpy as np

from matplotlib import pyplot as plt
plt.rcParams["animation.html"] = "jshtml" 
%matplotlib inline
from pprint import pprint
In [3]:
with open('./data/fish-preprocessed.pkl', 'rb') as f:
    data = pickle.load(f)
train_data = data['train']
test_data = data['test']

Zmienne kategoryczne w modelach liniowych

Podejść do "wektoryzacji" zmiennych kategorycznych jest wiele i ich stosowanie jest niezbędne w przypadku większości modeli. W przypadku modeli liniowych dobór odpowiedniego podejścia jest jednak wyjątkowo istotny, ponieważ model nie odtworzy ukrytych zależności - w przeciwieństwie do np. sieci neuronowych (przynajmniej w teorii). Oto wybór możliwości:

  1. Osobny model dla każdej z kategorii

    Mało wydajne, ponieważ zwiększamy liczbę stopni swobody (VC-dimensiom, klątwa wymiarowości itp) przy tej samej liczbie danych. Ale jeżeli każda z klas posiada inny proces generujący, to po co łączyć w jeden model?

  2. Utworzenie $k$ bądź $k - 1$ nowych zmiennych niezależnych (gdzie $k$ to liczba kategorii)

    Popularne i łatwe w użyciu podejście, nie wymaga zrozumienia zbioru danych, daje dobre wyniki dla modeli nie-liniowych. Można przeprowadzić na wiele sposobów - https://stats.idre.ucla.edu/spss/faq/coding-systems-for-categorical-variables-in-regression-analysis-2/

  3. Warunkowanie istniejących zmiennych niezależnych

    Wymaga zrozumienia zbioru danych (dla jakich zmiennych niezależnych wpływ na zmienną zależną może być różny zależnie od kategorii), ale dla modeli liniowych będzie dawać najlepsze wyniki (jeżeli przeprowadzone prawidłowo = zgodnie z procesem generatywnym/z jego najlepszą liniową aproksymacją).

    Poza trudnością w identyfikacji warunkowania, pojawiają się również trudności implementacyjne. Ręczne przygotowanie zbioru danych w ten sposób jest czasochłonne i łatwo o błąd. Bardzo dobrym rozwiązaniem jest podejście z pakietu R, tzw. formula function, którą posiada bardzo funkcjonalną składnię. W Pythonie z tej samej składni można skorzystać w pakiecie statsmodels.

Zastosujemy rozwiązanie nr 2. w wersji dummy / one-hot / one-of-K coding, czyli każda z klas będziała miała swoją kolumnę wypełnioną zerami i jedynkami. Dla jednej z klas (klasy referencyjnej) moglibyśmy nie tworzyć kolumny, ponieważ wyraz wolny "przejąłby" referencyjny wpływ na zmienną zależną. Gdybyśmy stosowali model inny niż liniowy z wyrazem wolnym takie podejście mogłoby się nie sprawdzić. Wyraz wolny, wkrtóce zniknie z naszego modelu liniowego, więc pozostawiamy $k$ kolumn.

In [4]:
one_hot = OneHotEncoder(drop=None, sparse=False)
train_data = train_data.join(pd.DataFrame(
    one_hot.fit_transform(train_data['Species'].to_numpy().reshape(-1, 1)),
    columns=[f'species_{i}' for i in range(7)],
    index=train_data.index,
))
train_data
Out[4]:
Species Weight Length3 Height Width intercept species_0 species_1 species_2 species_3 species_4 species_5 species_6
26 Bream 16.737549 40.6 16.3618 6.0900 1 1.0 0.0 0.0 0.0 0.0 0.0 0.0
141 Pike 19.806797 59.7 10.6863 6.9849 1 0.0 0.0 0.0 1.0 0.0 0.0 0.0
11 Bream 14.929327 36.2 14.3714 4.8146 1 1.0 0.0 0.0 0.0 0.0 0.0 0.0
69 Parkki 11.051308 25.8 10.3458 3.6636 1 0.0 1.0 0.0 0.0 0.0 0.0 0.0
23 Bream 16.443089 40.6 15.4686 6.1306 1 1.0 0.0 0.0 0.0 0.0 0.0 0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ...
72 Perch 2.234622 8.8 2.1120 1.4080 1 0.0 0.0 1.0 0.0 0.0 0.0 0.0
66 Parkki 9.766412 23.2 8.5376 3.2944 1 0.0 1.0 0.0 0.0 0.0 0.0 0.0
60 Whitefish 18.514038 43.5 12.3540 6.5250 1 0.0 0.0 0.0 0.0 0.0 0.0 1.0
94 Perch 10.006090 24.5 5.2185 3.6260 1 0.0 0.0 1.0 0.0 0.0 0.0 0.0
78 Perch 7.892379 19.4 5.1992 3.1234 1 0.0 0.0 1.0 0.0 0.0 0.0 0.0

110 rows × 13 columns

A teraz, wcześnie dopasowany na danych treningowych transformator należy zaaplikować na danych testowych:

In [5]:
test_data = test_data.join(pd.DataFrame(
    one_hot.transform(test_data['Species'].to_numpy().reshape(-1, 1)),
    columns=[f'species_{i}' for i in range(7)],
    index=test_data.index,
))

Oryginalną kolumnę zostawiamy w DataFrameach, ponieważ ułatwi nam to wizualizację.

Ridge Regression

Model: $$ y_i = \pmb{x}_i \pmb{\beta} $$

Regularyzacja $L_2$:

$$ \mathcal{R}(\pmb{\beta}) = \frac{\lambda}{2} \left\|\pmb{\beta}\right\|_2^2,$$

gdzie $\lambda$ to parametr metody regulujący siłę kary za złożoność.

Funkcja kosztu: $$ S(\pmb{\beta}) = \frac{1}{2} \sum_{i=1}^N \left| y_i - (\beta_0 + \pmb{x}_i \pmb{\beta}) \right|^2 + \frac{\lambda}{2} \left\|\pmb{\beta}\right\|^2, $$ gdzie macierz $X$ nie posiada tutaj dodatkowej kolumny $1$ odpowiadającej za wyraz wolny. Jest on traktowany osobno, ponieważ nie chcemy na niego zakładać regularyzacji.

Estymacja:

  1. Analityczna Najpierw musimy dopasować wyraz wolny $\beta_0$. Ma być to wartość o najmniejszym błędzie średniokwadratowym, gdyby wszystkie cechy miały wartość 0.

    $$ \hat{\beta}_0 = \frac{1}{N} \sum_{i=1}^N y_i $$

    Następnie wyraz wolny należy odjąć od zmiennej niezależnej co daje następującą funkcję kosztu: $$ S(\pmb{\beta}) = \frac{1}{2} \sum_{i=1}^N \left| (y_i - \hat{\beta}_0) - \pmb{x}_i \pmb{\beta}) \right|^2 + \frac{\lambda}{2} \left\|\pmb{\beta}\right\|^2, $$ którą można różniczkować względem parametrów $\pmb{\beta}$ i przyrównać do $0$. Dostajemy następujący estymator $$ \hat{\pmb{\beta}} = (\lambda I_N + X^\intercal X)^{-1}X^\intercal (Y - \hat{\beta}_0) $$

  2. Iteracyjna - Gradient Descent

    W tym podejściu do $S(\pmb{\beta})$ podchodzi się jak do funkcji kosztu nie wchodząc w jej interpretację i własności. Jest to problem z zakresu Convex Optimization.

Standaryzacja zmiennych

Standaryzowania zmiennych dokonuje się poprzez odjęcie średniej i przeskalowanie do jednostkowej wariancji.

Dla próbki $x$ ustandaryzowana wartość $z$ obliczana jest w następujący sposób:

$$ z = \frac{x - u}{s},$$

gdzie $u$ i $s$ to odpowiednio średnia i odchylenie standardowe populacji.

Dobrą praktyką jest standaryzować (ewentualnie skalować do przedziału) zmienne niezależne oraz zależne przed stosowaniem modeli. Dotychczas tego nie robiliśmy, ponieważ

  • dla modelu regresji liniowej skala zmiennych nie ma znaczenia
  • mieliśmy jeden model.

W poprzednim zeszycie wspomniane zostało, że przy regularyzowanej regresji skalowanie jest niezbędne - skoro nakładamy taką samą karę na wszystkie parametry, to powinny one mieć wpływ w tej same jednostce na wszystkie zmienne niezależne.

Ponadto dla ustandaryzowanych zmiennych łatwiej jest porównywać modele między sobą, a także dobierać rozkłady prior dla parametrów, a tym będziemy się zajmować w tym zeszycie.

Należy jednak pamiętać, że współczynniki dla ustandaryzowanych zmiennych są cięższe do interpretacji. Dlatego wiele pakietów statystycznych domyślnie dokonuje standaryzacji i raportuje parametry dla ustandaryzowanych zmiennych, a także przeskalowane do oryginalnych rozkładów cech.

In [6]:
train_data.columns
Out[6]:
Index(['Species', 'Weight', 'Length3', 'Height', 'Width', 'intercept',
       'species_0', 'species_1', 'species_2', 'species_3', 'species_4',
       'species_5', 'species_6'],
      dtype='object')

W tym zeszycie intercept nie będzie traktowany, jak kolejna z cech. Nie chcemy go także standaryzować - otrzymalibyśmy wektor zer.

Pewien niepokój powinno budzi standaryzowanie zmiennych powstałych w ramach one-hotowania. Pomimo, że mają wartości liczbowe, to ciężko o nich myśleć jako o zmiennych ciągłych. Z drugiej strony chcemy zachować tą samą skalę dla wszystkich zmiennych w celu regularyzacji. Nie istnieje jednoznacza odpowiedź na to pytanie. Większość debaty w Internecie sprowadza się do

sklearn zwróci błąd jeżeli wywołamy StandardScaler przed OneHotEncoder

Należy pamiętać, że jesteśmy teraz poza światem założeń, więc najlepiej zewaluować oba podejścia w różnych warunkach i na tej podstawie wyciągnąć wnioski.

In [7]:
features = [
    'Length3',
    'Height',
    'Width',
    'species_0',
    'species_1',
    'species_2',
    'species_3',
    'species_4',
    'species_5',
    'species_6',
]

target = 'Weight'
In [8]:
scaler = StandardScaler()
train_data[features] = scaler.fit_transform(train_data[features])
test_data[features] = scaler.transform(test_data[features])

Zastosujemy implementację Ridge Regression z pakietu sklearn, która estymuje wyraz wolny w pkt 1. powyżej, czyli nie regularyzuje go. Posiada ona następujące solvery (dokumentacja):

  • ‘auto’ chooses the solver automatically based on the type of data. ‘svd’ uses a Singular Value Decomposition of X to compute the Ridge coefficients. More stable for singular matrices than ‘cholesky’.

  • ‘cholesky’ uses the standard scipy.linalg.solve function to obtain a closed-form solution.

  • ‘sparse_cg’ uses the conjugate gradient solver as found in scipy.sparse.linalg.cg. As an iterative algorithm, this solver is more appropriate than ‘cholesky’ for large-scale data (possibility to set tol and max_iter).

  • ‘lsqr’ uses the dedicated regularized least-squares routine scipy.sparse.linalg.lsqr. It is the fastest and uses an iterative procedure.

  • ‘sag’ uses a Stochastic Average Gradient descent, and ‘saga’ uses its improved, unbiased version named SAGA. Both methods also use an iterative procedure, and are often faster than other solvers when both n_samples and n_features are large. Note that ‘sag’ and ‘saga’ fast convergence is only guaranteed on features with approximately the same scale. You can preprocess the data with a scaler from sklearn.preprocessing.

Wśród nich znajdziemy "analityczne" rozwiązanie w oparciu o operacje algebraiczne, a także interacyjne.

$\lambda$ nazywa się tutaj alpha :)

In [9]:
ridge_regression = linear_model.Ridge(
    alpha=1.0,
    fit_intercept=True,
    normalize=False
)

ridge_regression.fit(train_data[features], train_data[target])
print(ridge_regression.intercept_) # nie bedzie sie zmieniac wraz ze zmiana alpha
pprint(dict(zip(features, ridge_regression.coef_))) # gdy damy nieracjonalnie duze alpha to beda same zera
12.245810553153966
{'Height': 1.1263782709677466,
 'Length3': 2.5508051407598065,
 'Width': 1.2622090404664421,
 'species_0': -0.18340586231485526,
 'species_1': 0.08239548098583506,
 'species_2': 0.2285228104549227,
 'species_3': -0.07866971523456155,
 'species_4': 0.06262764460616334,
 'species_5': -0.2576554358615228,
 'species_6': 0.12506067436832283}
In [10]:
ridge_regression.score(test_data[features], test_data[target])
Out[10]:
0.9907396268991248

Od regularyzacji do Bayesian Linear Regression

(Zakładamy tutaj, że mamy już estymator wyrazu wolnego $\hat{\beta}_0$)

Przypomnijmy, że regresja liniowa to $$ (y_i - \hat{\beta}_0) = \pmb{x}_i \pmb{\beta} + \epsilon_i, $$ gdzie zakładamy $\epsilon_i \sim \mathcal{N}(0, \sigma^2)$. Czyli równoważnie możemy zapisać ten model jako: $$ (y_i - \hat{\beta}_0) | \pmb{x_i} \sim \mathcal{N}(\pmb{x}_i \pmb{\beta}, \sigma^2). $$

Log-likelihood takiego modelu to:

$$ \log p(D|\pmb{\beta}) = \frac{1}{N} \sum_{i=1}^N \log \left[ \frac{1}{\sqrt{2 \pi}\sigma}\exp \left( - \frac{((y_i - \hat{\beta}_0) - \pmb{x}_i \pmb{\beta})^2}{2 \sigma^2} \right)\right] = \text{const} - \frac{1}{2N\sigma^2}\sum_{i=1}^N \left| (y_i - \hat{\beta}_0) - \pmb{x}_i \pmb{\beta}) \right|^2 $$

Gdy spojrzymy na Ridge Regression od strony MAP (nawiązanie do ostatnich zajęć) $$ \arg\max_{\pmb{\beta}} \log p(\pmb{\beta}| D) = \arg\max_{\pmb{\beta}} [ \log p(\pmb{\beta}) + \log p (D | \pmb{\beta})], $$ to okazuje się, że Ridge Regression to probabilistyczny model ze standardowym Gaussowskim priorem na wagach.

Niech $\pmb{\beta} \sim \mathcal{N}(\pmb{m}, V)$: $$ \log p(\pmb{\beta}) = - \frac{1}{2} (\pmb{\beta} - \pmb{m})^\intercal V^{-1} (\pmb{\beta} - \pmb{m}) + \text{constant},$$ gdy $\pmb{m}=0$ i $V = I_{m}$ to $$ \log p(\pmb{\beta}) = - \frac{1}{2} \|\pmb{\beta}\|_2^2 + \text{constant},$$ czyli otrzymaliśmy Ridge Regression!

Bayesian Linear Regression

Mamy zatem taki sam model, jak w przypadku Linear Regression, ale dla parametrów wybieramy rozkład prior. Na początek załóżmy, że wariancja błędu jest znana i skupmy się współczynnikach. Likelihood modelu wyrażony jest w następujący sposób $$ p(Y|X, \pmb{\beta}, \beta_0, \sigma^2) = \mathcal{N}(Y | \beta_0 + X \pmb{\beta}, \sigma^2 \mathrm{I}_N) \propto \exp \left(- \frac{1}{2\sigma^2} (Y - \beta_0 \pmb{1}_N - X \pmb{\beta})^\intercal (Y - \beta_0 \pmb{1}_N - X \pmb{\beta}) \right) ,\tag{1}$$ gdzie przez $\mathcal{N}(Y | \beta_0 + X \pmb{\beta}, \sigma^2 \mathrm{I}_N)$ oznaczamy prawdopobieństwem obserwacji $Y$ przy rozkładzie normalnym o parametrach $(\beta_0 + X \pmb{\beta}, \sigma^2 \mathrm{I}_N)$.

Na parametr $\beta_0$ (wyraz wolny) zakładamy nieinformatywny prior, czyli jego estymacja sprowadza się do MLE i tak jak wcześniej $$ \hat{\beta}_0 = \frac{1}{N} \sum_{i=1}^N y_i = \bar{Y}.$$ Dla ułatwienia notacji będziemy zakładać, że $Y$ zostało wycentrowane (wtedy możemy pominąć wyraz wolny), zatem $$ Y \leftarrow Y - \bar{Y}.$$

Wyśrodkujmy zmienną niezależną w naszym zbiorze danych.

In [11]:
y_mean = train_data[target].mean()
train_data[target] = train_data[target] - y_mean
test_data[target] = test_data[target] - y_mean

Zatem $(1)$ upraszacza się do $$ p(Y|X, \pmb{\beta}, \sigma^2) \propto \exp \left(- \frac{1}{2\sigma^2} (Y - X \pmb{\beta})^\intercal (Y - X \pmb{\beta}) \right)$$ które w dalszym ciągu jest Gaussianem, więc rozkład sprzężony to również Gaussian, oznaczmy go w następujący sposób $$ p(\pmb{\beta}) = \mathcal{N}(\pmb{\beta} | \pmb{\beta}_0, \pmb{V}_0)$$

Wtedy posterior obliczamy jako $$ p(\pmb{\beta}| X, Y, \sigma^2) = \mathcal{N}(\pmb{\beta} | \pmb{\beta}_0, \pmb{V}_0) \mathcal{N}(Y | X\pmb{\beta}_0, \sigma^2 \mathrm{I}_N) = \mathcal{N}(\pmb{\beta} | \pmb{\beta}_N, \pmb{V}_N),\tag{2}$$ gdzie $$ \pmb{\beta}_N = \pmb{V}_N \pmb{V}_0^{-1} \pmb{\beta}_0+\frac{1}{\sigma^2}\pmb{V}_N X^\intercal Y $$ $$ \pmb{V}_N = \sigma^2 (\sigma^2 \pmb{V}_0^{-1} + X^\intercal X)^{-1} $$ $$ \pmb{V}_N^{-1} = \pmb{V}_0^{-1} + \frac{1}{\sigma^2} X^\intercal X $$

Skąd to się wzięło? Murphy:eq4.125, Bishop:eq2.116

  • Jeżeli $X$ jest ustandaryzowany, to rozsądnym jest $\pmb{\beta}_0 = \pmb{0}$
  • Jeżeli na tym etapie chcemy przestać być "Bayesowscy", to możemy zrobić MAP - tutaj znamy pełen rozkład (bo Gaussian), więc bierzemy modę (która równa jest średniej dla tego rozkładu)

Ale my nie chcemy opuszczać "Bayesowskiej" drogi i chcemy doprowadzić do znalezienia rozkładu predykcyjnego (ang. predictive distribution). Możemy tak zrobić, ponieważ znamy pełen posterior, a nie tylko funkcję do niego wprost proporcjonalną.

Przez $\mathcal{D}=(X, Y)$ oznaczmy dotychczasowy zbiór treningowy, składający się z $N$ par. Chcemy znaleźć rozkład (predykcyjny) $y$ dla nowej próbki $\pmb{x}$. W ogólność chcemy wykonać następującą operację $$ p(y|\pmb{x}, \mathcal{D}) = \int p(y|\pmb{x}, \pmb{\beta}) p(\pmb{\beta}|\mathcal{D}) \mathrm{d} \pmb{\beta},$$ gdzie $p(\pmb{\beta}|\mathcal{D})$ to posterior (policzony na bazie priora i likelihoodu). W naszym wypadku jest dodatkowo sparametryzowany przez $\sigma^2$.

Zatem finalnie

$$ p(y|\pmb{x}, \mathcal{D}, \sigma^2) = \int \mathcal{N}(y | \pmb{x}\pmb{\beta}, \sigma^2) \mathcal{N}(\pmb{\beta} | \pmb{\beta}_N, \pmb{V}_N) \mathrm{d} \pmb{\beta} = \mathcal{N}\left(y | x\pmb{\beta}_N, \sigma^2_N(\pmb{x})\right),\tag{3}$$

gdzie $$\sigma^2_N(\pmb{x}) = \sigma^2 + \pmb{x}\pmb{V}_N\pmb{x} $$

Warto zwrócić uwagę na $\sigma^2_N(\pmb{x})$, w której założona wariancja błędu jest powiększona o dodatkowy czynnik, który dodaje wariancji od posteriora wektora wag.

W przypadku probabilistycznego modelu nie mówimy już o przedziałach ufności i przedziałach predykcyjnych, tylko o credible intervals wynikających wprost z rozkładu predykcyjnego.

UWAGA 1: Taki model jest idealny dla uczenia przyrostowego. Dobry opis Bishop:p.154-156

UWAGA 2: Możemy znaleźć posterior dla modelu z priorem na $\sigma^2$ -> Murphy:p.236 rozdział 7.6.3 Bayesian inference when $\sigma^2$ is unknown

Empirical Bayes / Evidence Procedure

W powyższych rozważaniach mamy dwa hiper-parametry:

  • $\sigma^2$ - wariancja błędu; częściej mówi się o parametrze precision $\lambda=1/\sigma^2$
  • $\pmb{V}_0$ - macierz kowariancji dla priora na współczynnikach; najczęściej zakłamy sferyczny prior, więc $\pmb{V}_0 = \alpha^{-1} \mathrm{I}_N$ skupimy się na kolejnym precision $\alpha$

Oznaczmy $\eta= (\alpha, \lambda)$.

Popularnym rozwiązaniem tego problemu byłby Randomized Grid K-CV Search, gdzie dla losowego podzbioru konfiguracji z siatki przeprowadzona byłaby K-cross walidacja.

Alternatywą jest EB. Aby być w pełni Bayesowscy powinniśmy nałożyć hyperpriors na $\eta$. Wtedy możemy obliczyć posterior $\eta$ i rozkład predykcyjny ma następującą postać: $$ p(y|\pmb{x}, \mathcal{D}) = \int p(y|\pmb{x}, \pmb{\beta}) p(\pmb{\beta}|\mathcal{D}, \alpha, \lambda) p(\alpha, \lambda|\mathcal{D}) \mathrm{d} \pmb{\beta} \mathrm{d} \alpha \mathrm{d} \lambda$$

(Hiper-priory też mogą mieć swoje parametry i tak dalej...) Wycałkowanie $\alpha, \lambda$ jest analitycznie niemożliwe, obliczeniowo również. Pierwsze rozwiązanie to MAP - jeżeli posterior jest wyraźnie skupiony wokół wartości $\hat{\eta}$, to wybieramy ją i stosujemy jako plug-in estimate. By zaproponować alternatywę przypomnijmy jaką postać ma posterior $$ p(\alpha, \lambda | \mathcal{D}) \propto p(\mathcal{D}|\alpha, \lambda) p(\alpha, \lambda) $$ Jeżeli prior jest "płaski" i nie daje posteriora o skupionej masie, możemy zwrócić się ku MLE na likelihoodzie $p(\mathcal{D}|\alpha, \lambda)$. W obu przypadkach należy obliczyć posterior $$p(\mathcal{D}|\alpha, \lambda) = \int p(\mathcal{D}|\pmb{\beta}, \lambda) p(\pmb{\beta}|\alpha) \mathrm{d} \pmb{\beta},$$ a w przypadku MLE także przeprowadzić jego maksymalizację (opiera się na Expectation Maximization, które omówione zostanie w dalszej części semestru). Obie te rzeczy są dobrze opisane w Bishop:p.166-169 (UWAGA - inna notacja), rozdziały:

  • 3.5.1 Evaluation of the evidence function
  • 3.5.2 Maximizing the evidence function

UWAGA: W tym przypadku online learning nie jest wprost aplikowalny - "empiryczna część" potrzebuje próbki by dokonać maksymalizacji (posteriora bądź evidence). Można spróbować Expectation Maximization w konfiguracji batchowej - ciekawy temat na projekt.

In [12]:
from matplotlib.animation import FuncAnimation

class UpdateBayesianLinearRegression:
    '''Pomocnicza klasa do inkrementalnego uczenia'''
    def __init__(
        self,
        ax,
        model,
        X_test_with_cat,
        y_test,
        x_idx,
        x_label,
        y_label,
        category_idx=None,
        interval=0.95
    ):
        self.ax = ax
        self.model = model
        self.x_idx = x_idx
        
        self.interval = interval
        
        self.mean, = ax.plot([], [], color="red", label="Mean output")
        self.credible_interval = ax.fill_between(
            [],
            [],
            [],
            color='cornflowerblue',
            alpha=0.5,
            label=f"Credible Interval"
        )
        self.num_samples = ax.text(0.5, 1, '', fontsize=15, horizontalalignment='right',
                                   verticalalignment='bottom', transform=ax.transAxes)
        self.ratio_in_interval = ax.text(1, 0, '', fontsize=15, horizontalalignment='right',
                           verticalalignment='bottom', transform=ax.transAxes)
        
        if category_idx is not None:
            for cat in np.unique(X_test_with_cat[:, category_idx]):
                mask = X_test_with_cat[:, category_idx] == cat
                self.ax.scatter(
                    X_test_with_cat[mask, x_idx], y_test[mask],
                    marker='+', s=100,
                    alpha=1,
                    label=cat
                )
        else:
            self.ax.scatter(
                X_test_with_cat[:, x_idx], y_test,
                marker='+', s=100,
                alpha=1,
            )
        
        # Set up plot parameters
        self.ax.set_xlim(
            np.min(X_test_with_cat[:, x_idx])-0.2,
            np.max(X_test_with_cat[:, x_idx])+0.2
        )
        self.ax.set_xlabel(x_label)
        self.ax.set_ylabel(y_label)
        
        # Set up legend
        handles, labels = self.ax.get_legend_handles_labels()
        fig.legend(handles[:2], labels[:2], loc='upper right')
        if category_idx is not None:
            fig.legend(handles[2:], labels[2:], loc='lower right', title='Category')
        
    
    def model_fit(self, i, X_train, y_train):
        self.model.fit(X_train[i-1, :], y_train[i-1])

    def __call__(self, i, X_train, y_train, X_test, y_test):
        
        # This way the plot can continuously run and we just keep
        # watching new realizations of the process
        if i != 0:
            self.model_fit(i, X_train, y_train)

        predictive = self.model.predict(X_test)
        predictive_mean = predictive.mean()
        lower_interval, upper_interval = predictive.interval(self.interval)
        ci_ratio = np.sum(
            (lower_interval <= y_test) &
            (y_test <= upper_interval)
        ) / len(y_test)
        
        self.ax.set_ylim(
            min(np.min(lower_interval), np.min(y_test))-0.5,
            max(np.max(upper_interval), np.max(y_test))+0.5,
        )
        
        
        self.credible_interval.remove()
        self.credible_interval = self.ax.fill_between(
            X_test[:, self.x_idx],
            lower_interval,
            upper_interval,
            color='cornflowerblue',
            alpha=0.5,
            label=f"Credible Interval"
        )
        
        self.mean.set_data(X_test[:, self.x_idx], predictive_mean)
        self.num_samples.set_text(f'Number of train samples: {i}')
        self.ratio_in_interval.set_text(f'Percent of observations within CI: {ci_ratio*100:.2f}%')
        return self.mean, self.credible_interval, self.num_samples

Implementacja z rozkładem sprzężonym na $\pmb{\beta}$ i znanym $\alpha$ oraz $\lambda$

Zadanie 1 (2pkt)

W klasie poniżej zaimplementuj metodę predict tak by zwracała rozkład predykcyjny $Y$ dla zadanych danych wejściowych x. Wszystkie potrzebne wzory to $(2)$ i $(3)$, natomiast wszystkie potrzebne, nieznane parametry modelu zostały przypisane do zmiennych w konstruktorze oraz metodzie fit.

UWAGA: dodaj wyraz wolny do średniej

In [13]:
Student's answer Score: 2.0 / 2.0 (Top)
from scipy import stats

class ConjugateBayesLinReg:

    def __init__(self, n_features, alpha, lmbda, intercept=0):
        self.n_features = n_features
        self.alpha = alpha
        self.lmbda = lmbda
        self.mean = np.zeros(n_features)
        self.cov_inv = np.identity(n_features) * alpha
        self.cov = np.linalg.inv(self.cov_inv)
        self.intercept = intercept

    def fit(self, x, y):
        
        # If x and y are singletons, then we coerce them to a batch of length 1
        x = np.atleast_2d(x)
        y = np.atleast_1d(y)

        # Update the inverse covariance matrix (Bishop eq. 3.51)
        cov_inv = self.cov_inv + self.lmbda * x.T @ x # Vn-1

        # Update the mean vector (Bishop eq. 3.50)
        cov = np.linalg.inv(cov_inv) # Vn
        mean = cov @ (self.cov_inv @ self.mean + self.lmbda * y @ x) #beta_n

        self.cov_inv = cov_inv
        self.cov = cov
        self.mean = mean


    def predict(self, x):
        '''Metoda zwraca rozkład predykcyjny zmiennej niezależnej na bazie bieżących parametrów modelu.'''
        # If x and y are singletons, then we coerce them to a batch of length 1
        x = np.atleast_2d(x)
        
        # TU WPISZ KOD !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        sigma = 1 / self.lmbda
        y_pred_var = np.diagonal(sigma + x @ self.cov @ x.T)# sigma^2 + x Vn x
        y_pred_mean = x @ self.mean + self.intercept 

        # Drop a dimension from the mean and variance in case x and y were singletons
        y_pred_mean = np.squeeze(y_pred_mean)
        y_pred_var = np.squeeze(y_pred_var)
    
        return stats.norm(loc=y_pred_mean, scale=y_pred_var ** .5)

    @property
    def weights_dist(self):
        return stats.multivariate_normal(mean=self.mean, cov=self.cov)
In [14]:
fig, ax = plt.subplots(figsize=(16, 8))
ud = UpdateBayesianLinearRegression(
    ax,
    model=ConjugateBayesLinReg(n_features=len(features), alpha=.3, lmbda=1),
    X_test_with_cat=test_data[features+['Species']].to_numpy(),
    x_idx=1,
    x_label="Length3",
    y_label="Weight",
    category_idx=-1,
    y_test=test_data[target],
)

anim = FuncAnimation(
    fig, 
    ud,
    fargs=(
        train_data[features].to_numpy(),
        train_data[target].to_numpy(),
        test_data[features].sort_values(by=features[ud.x_idx], axis=0).to_numpy(),
        test_data[[target, features[ud.x_idx]]].sort_values(by=features[ud.x_idx], axis=0)[target].to_numpy()
    ),
    frames=train_data.shape[0],
    interval=100,
    blit=True,
)
plt.close()
anim
Out[14]:
In [15]:
pprint(dict(zip(features, ud.model.mean)))
pprint(ud.model.cov)
{'Height': 1.0862281243634957,
 'Length3': 2.796021313804139,
 'Width': 1.109668534622248,
 'species_0': -0.17910669771856647,
 'species_1': 0.1019449027684658,
 'species_2': 0.26136563778587174,
 'species_3': -0.18130359991279477,
 'species_4': 0.07624150094483184,
 'species_5': -0.246153085303078,
 'species_6': 0.13430952784430872}
array([[ 0.24667031, -0.13081177, -0.10277793,  0.04713826,  0.02911855,
         0.01412368, -0.11700408,  0.00614219, -0.00217489,  0.00711556],
       [-0.13081177,  0.43262431, -0.19006008, -0.1805213 , -0.06013925,
         0.06120188,  0.12771946,  0.02174284,  0.04344136,  0.00840203],
       [-0.10277793, -0.19006008,  0.22520651,  0.08077008,  0.02079129,
        -0.05449351,  0.00175367, -0.01673434, -0.01852415, -0.01340347],
       [ 0.04713826, -0.1805213 ,  0.08077008,  0.80568831,  0.47778094,
         0.81035292,  0.49499441,  0.55398073,  0.48163996,  0.32374973],
       [ 0.02911855, -0.06013925,  0.02079129,  0.47778094,  0.30299822,
         0.51858843,  0.32007709,  0.3526548 ,  0.31172131,  0.20454694],
       [ 0.01412368,  0.06120188, -0.05449351,  0.81035292,  0.51858843,
         0.98931378,  0.63731075,  0.65837643,  0.58871251,  0.38195204],
       [-0.11700408,  0.12771946,  0.00175367,  0.49499441,  0.32007709,
         0.63731075,  0.48699352,  0.42690305,  0.3863347 ,  0.2459788 ],
       [ 0.00614219,  0.02174284, -0.01673434,  0.55398073,  0.3526548 ,
         0.65837643,  0.42690305,  0.45111044,  0.39675893,  0.25625006],
       [-0.00217489,  0.04344136, -0.01852415,  0.48163996,  0.31172131,
         0.58871251,  0.3863347 ,  0.39675893,  0.36734145,  0.227785  ],
       [ 0.00711556,  0.00840203, -0.01340347,  0.32374973,  0.20454694,
         0.38195204,  0.2459788 ,  0.25625006,  0.227785  ,  0.15779719]])

Implementacja z rozkładem sprzężonym na $\pmb{\beta}$ i nieznanym $\alpha$ oraz $\lambda$

Jak wspomnieliśmy powyżej ta metoda może opierać się na:

  1. MLE na likelihoodzie (Evidence Maximization) Wyprowadzenie (na podstawie Bishop) oraz implementację tego podejścia można znaleźć pod tym linkiem - http://krasserm.github.io/2019/02/23/bayesian-linear-regression/
  2. MAP Rozkładem sprzężonym dla precision jest rozkład Gamma, który zależnie od parametrów preferuje małe bądź większe wartości. Implementacja tego podejścia jest dostępna w pakiecie sklearn pod nazwą sklearn.linear_model.BayesianRidge. Wykorzystamy ją tutaj wraz z omówieniem parametryzacji.

    UWAGA: tutaj $\lambda$ i $\alpha$ są zamienione!

In [16]:
from sklearn.linear_model import BayesianRidge

Wybrane parametry klasy BayesianRidge:

  • n_iter: int, default=300

    Maximum number of iterations. Should be greater than or equal to 1.

  • tol: float, default=1e-3

    Stop the algorithm if w has converged.

  • alpha_1: float, default=1e-6

    Hyper-parameter : shape parameter for the Gamma distribution prior over the alpha parameter.

  • alpha_2: float, default=1e-6

    Hyper-parameter : inverse scale parameter (rate parameter) for the Gamma distribution prior over the alpha parameter.

  • lambda_1: float, default=1e-6

    Hyper-parameter : shape parameter for the Gamma distribution prior over the lambda parameter.

  • lambda_2: float, default=1e-6

    Hyper-parameter : inverse scale parameter (rate parameter) for the Gamma distribution prior over the lambda parameter.

  • `alpha_init: float, default=None

    Initial value for alpha (precision of the noise). If not set, alpha_init is 1/Var(y).

       New in version 0.22.
  • lambda_init: float, default=None

    Initial value for lambda (precision of the weights). If not set, lambda_init is 1.

       New in version 0.22.
  • fit_intercept: bool, default=True

    Whether to calculate the intercept for this model. The intercept is not treated as a probabilistic parameter and thus has no associated variance. If set to False, no intercept will be used in calculations (i.e. data is expected to be centered).

Ta implementacja nie wspiera inkrementalnego przetwarzania (w sklearn metody działające w trybie online posiadają metodę partial_fit), więc w każdym kroku będziemy wykorzystywać narastające okno danych. Ponadto jako inicjalną wartość $\lambda$ i $\alpha$ będziemy podawać estymatę z poprzedniego okna.

In [17]:
class UpdateBayesianLinearRegressionInreasingWindow(UpdateBayesianLinearRegression):
    def model_fit(self, i, X_train, y_train):
        self.model.fit(X_train[:i, :], y_train[:i])

Zdefiniujemy teraz klasę FullConjugateBayesLinReg na bazie implementacji z sklearn - sprametryzujemy ją zgodnie z naszą notacją.

In [18]:
from scipy import stats

class FullConjugateBayesLinReg(ConjugateBayesLinReg):

    def __init__(
        self,
        n_features,
        n_iter=300,
        tol=1.e-3,
        alpha_shape=1.e-6, # non-informative prior 
        alpha_rate=1.e-6,
        lmbda_shape=2, # mode=1 => standard normal posterior
        lmbda_rate=1,
        alpha_init=None,
        lmbda_init=1,
        fit_intercept=False
    ):
        self.n_iter = n_iter
        self.tol = tol
        self.alpha_shape = alpha_shape
        self.alpha_rate = alpha_rate
        self.lmbda_shape = lmbda_shape
        self.lmbda_rate = lmbda_rate
        self.alpha = alpha_init
        self.lmbda = lmbda_init
        self.fit_intercept = fit_intercept
        
        self.mean = np.zeros((n_features,))
        self.cov = np.eye(n_features) / self.alpha
        self.intercept = 0

    def fit(self, x, y):
        
        # If x and y are singletons, then we coerce them to a batch of length 1
        x = np.atleast_2d(x)
        y = np.atleast_1d(y)

        model = BayesianRidge(
            n_iter=self.n_iter,
            tol=self.tol,
            alpha_1=self.lmbda_shape,
            alpha_2=self.lmbda_rate,
            lambda_1=self.alpha_shape,
            lambda_2=self.alpha_rate,
            alpha_init=self.lmbda,
            lambda_init=self.alpha,
            fit_intercept=self.fit_intercept,
            normalize=False,
        )
        
        model.fit(x, y)
        
        self.mean = model.coef_
        self.cov = model.sigma_
        self.alpha = model.lambda_
        self.lmbda = model.alpha_
        self.intercept = model.intercept_
In [19]:
fig, ax = plt.subplots(figsize=(16, 8))
ud = UpdateBayesianLinearRegressionInreasingWindow(
    ax,
    model=FullConjugateBayesLinReg(n_features=len(features), alpha_init=.3, lmbda_init=1),
    X_test_with_cat=test_data[features+['Species']].to_numpy(),
    x_idx=1,
    x_label="Length3",
    y_label="Weight",
    category_idx=-1,
    y_test=test_data[target],
)

anim = FuncAnimation(
    fig, 
    ud,
    fargs=(
        train_data[features].to_numpy(),
        train_data[target].to_numpy(),
        test_data[features].sort_values(by=features[ud.x_idx], axis=0).to_numpy(),
        test_data[[target, features[ud.x_idx]]].sort_values(by=features[ud.x_idx], axis=0)[target].to_numpy()
    ),
    frames=train_data.shape[0],
    interval=100,
    blit=True,
)
plt.close()
anim
Out[19]:
In [20]:
pprint(dict(zip(features, ud.model.mean)))
print('lambda:', ud.model.lmbda)
print('alpha:', ud.model.alpha)
pprint(ud.model.cov)
{'Height': 1.0706431607315405,
 'Length3': 2.865785514978405,
 'Width': 1.068739728400549,
 'species_0': -0.17572775204966468,
 'species_1': 0.10809381994597501,
 'species_2': 0.2693189978625679,
 'species_3': -0.21116165465397435,
 'species_4': 0.07993841676750216,
 'species_5': -0.243253026378295,
 'species_6': 0.13689589044451547}
lambda: 5.2761397393466885
alpha: 0.8197761008230083
array([[ 4.98574760e-02, -2.76479556e-02, -2.00075200e-02,
         1.00671033e-02,  6.03423923e-03,  2.63537432e-03,
        -2.38876411e-02,  1.15295811e-03, -6.04678963e-04,
         1.40566864e-03],
       [-2.76479556e-02,  9.11102453e-02, -4.00361544e-02,
        -3.80322789e-02, -1.26998432e-02,  1.28959605e-02,
         2.69800727e-02,  4.56514032e-03,  9.11015229e-03,
         1.77793249e-03],
       [-2.00075200e-02, -4.00361544e-02,  4.59016680e-02,
         1.70370897e-02,  4.42893407e-03, -1.11947789e-02,
        -2.10797198e-04, -3.46852288e-03, -3.95636434e-03,
        -2.71665166e-03],
       [ 1.00671033e-02, -3.80322789e-02,  1.70370897e-02,
         2.81938816e-01,  1.71346781e-01,  3.00866044e-01,
         1.89020697e-01,  2.04567779e-01,  1.79753324e-01,
         1.19128250e-01],
       [ 6.03423923e-03, -1.26998432e-02,  4.42893407e-03,
         1.71346781e-01,  1.08080260e-01,  1.91072748e-01,
         1.20802893e-01,  1.29509069e-01,  1.14820399e-01,
         7.51197169e-02],
       [ 2.63537432e-03,  1.28959605e-02, -1.11947789e-02,
         3.00866044e-01,  1.91072748e-01,  3.58979343e-01,
         2.32633731e-01,  2.40416847e-01,  2.14615879e-01,
         1.39440564e-01],
       [-2.38876411e-02,  2.69800727e-02, -2.10797198e-04,
         1.89020697e-01,  1.20802893e-01,  2.32633731e-01,
         1.66188382e-01,  1.56300240e-01,  1.40514838e-01,
         9.02989156e-02],
       [ 1.15295811e-03,  4.56514032e-03, -3.46852288e-03,
         2.04567779e-01,  1.29509069e-01,  2.40416847e-01,
         1.56300240e-01,  1.63507541e-01,  1.44686753e-01,
         9.37985890e-02],
       [-6.04678963e-04,  9.11015229e-03, -3.95636434e-03,
         1.79753324e-01,  1.14820399e-01,  2.14615879e-01,
         1.40514838e-01,  1.44686753e-01,  1.31553802e-01,
         8.34588535e-02],
       [ 1.40566864e-03,  1.77793249e-03, -2.71665166e-03,
         1.19128250e-01,  7.51197169e-02,  1.39440564e-01,
         9.02989156e-02,  9.37985890e-02,  8.34588535e-02,
         5.61386065e-02]])

Zadanie 2

Pewien student posiada samochód i dosyć często podróżuje z Wrocławia w Bieszczady. Jest oszczędny, więc wybiera drogę z pominięciem odcinków płatnych. Student (ani my) nie wie ile kilometrów ma ta droga, nie zna spalania swojego samochodu, a chciałby zaplanować budżet na pewną wyprawę. W ciągu ostatnich 10 podróży zanotował następujące dane:

  • pieniądze wydane na paliwo [PLN] - zkładamy, że we Wrocławiu jest na stacji z pustym bakiem i dojeżdża w góry na ostatnich oparach
  • średnia prędkość jazdy w danej podróży [km/h] - uwzględnia tylko czas przemieszczania się
  • liczba postojów podczas danej podróży - nie zanotował jednak ile postoje trwały
  • łączny czas podróży [h]

W całym rozumowaniu będziemy mieć jedno, istotne i naciągane założenie - spalanie samochodu nie zależy od prędkości jazdy i nie musi być stałe w każdej podróży.

Owy student podjął decyzję, że należy rzucić wszystko, zatankować samochód za 100 PLN i stając dwa razy na kawę jechać 90km/h aż do skończenia paliwa. O tej porze roku zachód słońca jest o godzinie 17:30. Ile maksymalnie godzin przed 17:30 powinien wyjechać żeby z 90% szansą zatrzymać się o tej porze, bądź póżniej?

Zastosuj metody poznane na zajęciach aby wyznaczyć szukaną wartość.

In [21]:
df = pd.read_pickle('./data/fuel.pkl')
df
Out[21]:
money speed num_stops time
0 245.622776 103.635372 3 5.445866
1 233.260688 78.237048 1 6.684284
2 231.230264 118.704676 2 5.239483
3 268.520137 101.427849 2 6.215113
4 246.789517 53.266307 4 11.231389
5 237.886947 94.934671 4 6.285680
6 246.226771 68.733270 5 9.192864
7 194.570072 85.653516 5 5.992770
8 230.535203 96.105184 4 7.258456
9 216.144688 69.252206 2 7.525162

Zadanie 2.1 (2 pkt)

W komórce poniżej przygotuj dane z df do zaaplikowania następującego modelu liniowego: \begin{equation} \text{time} = \frac{\text{money}}{\text{speed}} * \text{distance_for_money_unit} + \text{num_stops} * \text{stop_time}, \end{equation}

gdzie:

  • $\text{distance_for_money_unit}$ to współczynnik odpowiadający za dystans przejechany za daną kwotę
  • $\text{stop_time}$ to współczynnik odpowiadający za czas spędzony na postoju
  1. Utwórz zmienną 'money/speed' w zbiorze danych
  2. Dokonaj standaryzacji zmiennych niezależnych, które będą wykorzystane w modelu
  3. Wycentruj zmienną niezależną i zapisz jej średnią pod zmienną y_mean
In [22]:
Student's answer Score: 1.0 / 2.0 (Top)
features = [
    'money/speed',
    'num_stops'
]
target = 'time'
scaler = StandardScaler()
y_mean = df[target].mean()

# TU WPISZ KOD

df[features[0]] = df["money"] / df["speed"] 

df[features] = scaler.fit_transform(df[features])
df[target] = df[target] - y_mean
df.head()
Out[22]:
money speed num_stops time money/speed
0 245.622776 103.635372 -0.150756 -1.661241 -0.638906
1 233.260688 78.237048 -1.658312 -0.422823 0.181907
2 231.230264 118.704676 -0.904534 -1.867624 -1.205615
3 268.520137 101.427849 -0.904534 -0.891994 -0.266578
4 246.789517 53.266307 0.603023 4.124282 2.399314

Zadanie 2.2 (1pkt)

W komórce poniżej zainicjalizuj wybrany model Bayesowskiej Regresji Liniowej pod zmienną fuel_model i wykorzystaj go w animacji inkrementacyjnego uczenia - animacja jest w pełni przygotowana.

1.ConjugateBayesLinReg

In [23]:
Student's answer Score: 1.0 / 0.0 (Top)
# TU WPISZ KOD
fuel_model = ConjugateBayesLinReg(n_features=len(features), alpha=.3, lmbda=1)



fig, ax = plt.subplots(figsize=(16, 8))
ud = UpdateBayesianLinearRegressionInreasingWindow(
    ax,
    model=fuel_model,
    X_test_with_cat=df[features].to_numpy(),
    x_idx=0,
    x_label='money/speed',
    y_label="Travel time",
    y_test=df[target],
)

anim = FuncAnimation(
    fig, 
    ud,
    fargs=(
        df[features].to_numpy(),
        df[target].to_numpy(),
        df[features].sort_values(by=features[ud.x_idx], axis=0).to_numpy(),
        df[[target, features[ud.x_idx]]].sort_values(by=features[ud.x_idx], axis=0)[target].to_numpy()
    ),
    frames=df.shape[0],
    interval=100,
    blit=True,
)
plt.close()
anim
Out[23]:

2. FullConjugateBayesLinReg

In [24]:
# TU WPISZ KOD
fuel_model = FullConjugateBayesLinReg(n_features=len(features), alpha_init=.3, lmbda_init=1)



fig, ax = plt.subplots(figsize=(16, 8))
ud = UpdateBayesianLinearRegressionInreasingWindow(
    ax,
    model=fuel_model,
    X_test_with_cat=df[features].to_numpy(),
    x_idx=0,
    x_label='money/speed',
    y_label="Travel time",
    y_test=df[target],
)

anim = FuncAnimation(
    fig, 
    ud,
    fargs=(
        df[features].to_numpy(),
        df[target].to_numpy(),
        df[features].sort_values(by=features[ud.x_idx], axis=0).to_numpy(),
        df[[target, features[ud.x_idx]]].sort_values(by=features[ud.x_idx], axis=0)[target].to_numpy()
    ),
    frames=df.shape[0],
    interval=100,
    blit=True,
)
plt.close()
anim
Out[24]:

Dane testowe na podstawie treści zadania:

  • 100 PLN
  • 90km/h
  • 2 postoje
In [25]:
test_data = np.array([100/90, 2]).reshape(1, -1)

Zadanie 2.3 (1pkt)

Szukamy teraz wartości od, której $\text{time}$ jest większy z prawdopobieństwem co najmniej 0.9 (90% szans na dojechanie). Równoważnie wartości, od której $\text{time}$ jest mniejszy z prawdopobieństwem co najmniej 0.1. Czyli kwantylu rzędu $p=0.1$.

$$P_X((-\infty, x_p])) \geqslant p \iff P_X([x_p, \infty))) \geqslant 1-p$$

Dla rozkładów z modułu stats pakietu scipy ta funkcjonalność zaimplementowana jest pod metodą ppf.

Na podstawie rozkładu predykcyjnego dla danych testowych test_data oblicz maksymalną liczbę godzin (w sytemie dziesiętnym - 15min = 0.25h itd) zgodnie z treścią zadania i przypisz wartość do zmiennej worst_case_time.

UWAGA: pamiętaj jakie przekształcenia były stosowane na zbiorze danych!

In [26]:
Student's answer Score: 1.25 / 1.0 (Top)
# TU WPISZ KOD
test_data_scaled = scaler.transform(test_data)
y_pred = fuel_model.predict(test_data_scaled)

fig, ax = plt.subplots()
x= np.arange(-10,4,0.001)
ax.set_title(f'N({y_pred.mean():.2f},{y_pred.std():.2f})')
ax.set_xlabel('x')
ax.set_ylabel('f(x)')
ax.plot(x, y_pred.pdf(x))

res = y_pred.ppf(q=0.1)
px=np.arange(-10,res,0.01)
ax.fill_between(px,y_pred.pdf(px),alpha=0.5, color='g')

plt.show()

worst_case_time = res +  y_mean

hours = np.floor(worst_case_time)
minutes = (worst_case_time - hours) * 60
minutes_floor = np.floor(minutes)
seconds = round(minutes - minutes_floor,2) * 60
print(f'{hours:.0f}h {minutes_floor:.0f}m {seconds}s')
print(f'worst time  = {worst_case_time:.2f}')
2h 14m 25.2s
worst time  = 2.24

Implementacja z dowolonym rozkładem na $\pmb{\beta}$, $\alpha$ oraz $\lambda$ z wykorzystaniem Variational Inference

W sytuacji gdy chcemy operować na dowolnych rozkładach prior nie jesteśmy w stanie przeprowadzić dokładnej estymacji Bayesowskiej. Ewaluacja rozkładów posterior w sposób numeryczny jest w ogólności poza zasięgiem możliwości obliczeniowych. Rozwiązaniem tego problemu jest przybliżone wnioskowanie Bayesowskie (ang. approximate Bayesian inference), gdzie miedzy innymi możemy wymienić następujące metody:

  • Markov Chain Monte Carlo - poprzez błądzie losowe w oparciu o Markov Chain daje aproksymację skomplikowanych rozkładów ciągłych za pomocą techniki Monte Carlo
  • Expectation Maximization - w interacyjnej procedurze odszukuje lokalne optima parametrów dla MLE i MAP gdy model zależy od pewnych nieobserwowanych zmiennych
  • Variational Inference - problem maksymalizacji prawdopobieństwa brzegowego (Marginal Likelihood) jest przeniesiony na problem optymalizacji gradientowej na dolnym ograniczeniu likelihoodu - Evidence Lower BOund (ELBO)

Biblioteka pyro wspiera MCMC oraz VI, a konkretnie rozszerzoną wersję Stochastic Variational Inference, która umożliwia przetwrzania batchowe - analogicznie do Stochastic Gradient Descent.